Physics-informed attention-based neural network for hyperbolic partial differential equations: application to the Buckley–Leverett problem

Physics-informed neural networks (PINNs) have enabled significant improvements in modelling physical processes described by partial differential equations (PDEs) and are in principle capable of modeling a large variety of differential equations. PINNs are based on simple architectures, and learn the behavior of complex physical systems by optimizing the network parameters to minimize the residual of the underlying PDE. Current network architectures share some of the limitations of classical numerical discretization schemes when applied to non-linear differential equations in continuum mechanics. A paradigmatic example is the solution of hyperbolic conservation laws that develop highly localized nonlinear shock waves. Learning solutions of PDEs with dominant hyperbolic character is a challenge for current PINN approaches, which rely, like most grid-based numerical schemes, on adding artificial dissipation. Here, we address the fundamental question of which network architectures are best suited to learn the complex behavior of non-linear PDEs. We focus on network architecture rather than on residual regularization. Our new methodology, called physics-informed attention-based neural networks (PIANNs), is a combination of recurrent neural networks and attention mechanisms. The attention mechanism adapts the behavior of the deep neural network to the non-linear features of the solution, and break the current limitations of PINNs. We find that PIANNs effectively capture the shock front in a hyperbolic model problem, and are capable of providing high-quality solutions inside the convex hull of the training set.

www.nature.com/scientificreports/ be used to assimilate data and observations into numerical models, or be used in parameter identification (the inverse problem) 19,38 and uncertainty quantification [39][40][41][42] . Learning solutions of nonlinear PDEs using current network architectures presents some of the same limitations of classical numerical discretization schemes. For instance, one of the main limitations of classical numerical methods (e.g. finite differences, volumes, elements) is the need to devise suitable upwind discretizations that yield smooth and accurate solutions near the shock fronts. Methods that use centered approximations, without numerical dissipation, such as standard finite differences or the Galerkin finite element method, lead to solutions that are polluted by spurious oscillations around the shock front, often leading to numerical instabilities. PINN might present similar issues if the same schemes are used to calculate residuals.
A paradigmatic example is the solution of hyperbolic PDEs. Hyperbolic conservation laws describe a plethora of physical systems in gas dynamics, acoustics, elastodynamics, optics, geophysics, and biomechanics 43 . Hyperbolic PDEs are challenging to solve numerically using classical discretization schemes, because they tend to form self-sharpening, highly-localized, nonlinear shock waves that require specific approximation strategies and fine meshes 44 . Solving hyperbolic PDEs seems to be challenging for neural networks as well, as the ability of current PINNs to learn PDEs with a dominant hyperbolic character relies on adding artificial dissipation [45][46][47][48] , or on using a priori knowledge to increase the number of training points along the shock trajectories 39 or adaptive activation functions 49,50 .
In this work, we propose a new perspective on solving hyperbolic PDEs and traditional limitations of classical numerical methods using deep learning. The proposed method relies in two core ideas: (1) A modified PINN architecture can provide a more general method for solving hyperbolic conservation law problems without a priori knowledge or residual regularization. (2) Relating network architecture with the physics encapsulated in a given PDE is possible and has a beneficial impact. Our hypothesis is that sophisticated, physics-specific network architectures (i.e. networks whose internal hierarchy is related to the physical processes being learned) may be more effectively trained and easier understood than standard feed-forward multilayer perceptrons.
The proposed new architecture uses attention mechanisms to automatically detect shocks in the solution of hyperbolic PDEs. Our use of attention mechanisms to enrich PINN network architectures is inspired by recent advances in deep learning for language processing and translation 51,52 . The resulting network is a combination of gated recurrent units (GRUs) and attention mechanisms; we call this a physics-informed attention-based neural network (PIANN). The combination of both elements in the architecture allows for determination of the most relevant information (recurrent neural network with memory) to adapt the behavior of the deep neural network to approximate sharp shocks without the necessity of residual regularization or a priori knowledge (attention mechanism).
Previous works as 19,46,53 , introduced initial and boundary conditions in the formulation as a penalty term in the objective function. The main drawback of this approach is, if this term was not exactly zero after training, the boundary condition is not completely satisfied. Other authors as Lagaris et al. 54 or more recently, Schiassi et al. 55,56 make use of the Extreme Theory of Functional Connections 57 to enforce the initial and boundary conditions in the solution. As in these works, we also enforce the initial and boundary conditions as hard constraints.
We test our PIANN approach by solving a classical hyperbolic model problem, namely the Buckley-Leverett equation 44,58 . The Buckley-Leverett equation with non-convex flux function is an excellent benchmark to test the overall potential of PIANNs in solving hyperbolic PDEs. We find that PIANNs effectively capture the shock front propagation and are capable of providing high quality solutions for mobility ratios inside the convex hull of training set. Remarkably, PIANNs are able to provide smooth, accurate shock fronts without explicitly introducing additional constraints or dissipation in the residual, through an artificial diffusion term or upwinding the spatial derivatives.
The rest of the paper is organized as follows. The mathematical formulation of the problem, describing incompressible two-phase flow in porous media is introduced in "Section Problem formulation". "Section Methodology" introduces our methodology to address the problem using deep neural networks. The proposed PIANN architecture is described in "Section PIANN architecture". Experimental results are reported in "Section Results" and finally, "Section Discussion" concludes the paper.

Problem formulation
The problem of interest is that of two immiscible fluids flowing through a horizontal porous medium. We further assume that the fluids and overall system are incompressible. The Buckley-Leverett (BL) equation 58  This first-order hyperbolic equation is of interest as its solution can display both smooth solutions (rarefactions) and sharp fronts (shocks). Although the solution to this problem can be obtained analytically for simple (2) u M (x, 0) = 0, ∀x > 0, Initial condition www.nature.com/scientificreports/ one-dimensional settings, the precise and stable resolution of these shocks poses well-known challenges for numerical methods 44 . Physics-informed neural networks (PINNs) have been tested on this problem by Fuks and Tchelepi 45 who report good performance for concave fractional flow functions. In addition, Fraces and Tchelepi 59 provide an accurate solution introducing two physical constrains and subsequently modifying the fractional flux equation by a piece-wise form which differentiable form allow to capture the shock. However, these constraints are problem dependent and do not provide a general solution for hyperbolic equations in heterogeneous media. Whether the solution of the BL problem with non-convex flux function can be learnt by deep neural networks without the aid of artificial physical constraints remains an open question.
We take f M to be the S-shaped flux function for which we can obtain the analytical solution of the problem: , and u * represents the shock location defined by the Rankine-Hugoniot condition 44 . Note that equation (5) describes a family of functions of space, time, characterized by fluid mobility ratio M, u M (x, t) . In "Section Results" we use this analytical solution to test the accuracy of the PIANN model.

Methodology
T be a discrete version of the domain of u M . We define our PIANN as a vector function u θ : where θ are the weights of the network to be estimated during training. The inputs for the proposed architecture are pairs of (t, M) and the output is a vector where the i-th component is the solution evaluated in x i . Notice the different treatment applied to spatial and temporal coordinates. Whereas t is a variable of the vector function u θ , the locations where we calculated the solution x 0 , . . . , x N are fixed in advance. The output is a saturation map and therefore its values have to be in the interval [0, 1].
For the sake of readability, we introduce the architecture of u θ in "section PIANN architecture". However, we note already here that in order to enforce the boundary condition, we let our PIANN learn only the components u θ (t, M) 1 , . . After utilizing the information provided by the initial and boundary conditions enforcing u θ (0, M) i = 0, i = 1, . . . , N and u θ (t, M) 0 = 1 , respectively, we now define a loss function based on the information provided by Eq. (1). To calculate the first term we propose two options. The first option is a central finite difference approximation, that is, Alternatively, we can calculate the derivative of our PIANN with respect to t since we know the functional form of u θ . It can be calculated using the automatic differentiation tools included in many machine learning libraries, such as Pytorch. Thus, we propose a second option to calculate this term as The second term of Eq. (1), the derivative of the flux with respect to the spatial coordinate, is approximated using central finite difference as where the vector of fluxes at the i-th location x i is calculated as The spatial coordinate x is included as a fixed parameter in our architecture.
The loss function to estimate the parameters of the PINN is given as www.nature.com/scientificreports/ where � · � F is the Frobenius norm. There are two main lines in PINNs literature to deal with the initial and boundary conditions. The first one is to include them as a penalty term in the objective function 19,46,53 . The main drawback of this approach is, if this penalty term is not zero after training, the initial and boundary conditions are not totally satisfied by the solution.
To solve this problem, we directly enforce the initial and boundary conditions in the architecture of our PIANN. Thus, we follow the path of the second line in PINN literature [54][55][56] that enforces initial and boundary conditions as hard constraints. With respect to the first line, this provides several advantages. First, we enforce a stronger constraint that does not allow any error on the initial and boundary conditions. Second, the PIANN does not need to learn these conditions by itself, and it can instead concentrate exclusively on learning the parameters that minimize the residuals of the BL equation. Third, there are no weights to be tuned to control the effect of the initial and boundary conditions in the final solution.
The parameters of the PIANN are estimated using the Adam optimizer 60 to minimize Eq. (9) with respect to θ . Training algorithm pseudo-code is provided in Algorithm 1. Calculate the first term of the residual R 1 (T , M) as in eq. (6). 4: Calculate the values of the flux function f T (t, M) using eq. (8).

5:
Calculate the second term of the residual R 2 (T , M) using the values of flux function in eq. (7). 6: Calculate the objective function L (T ) as in eq. (9). 7: Calculate the derivatives of L (T ) with respect to T using automatic differentiation. 8: Update the parameters T of the PIANN according to Adam optimizer steps.

PIANN architecture
Although it has been demonstrated that neural networks are universal function approximators, the proper choice of architecture generally aids learning and certain challenging problems (e.g. solving non-linear PDEs) may require more specific architectures to capture all their properties (see for instance Extreme Learning Machine 61,62 ). For that reason, we have proposed a new architecture, inspired by 51 , to solve non-linear PDEs with discontinuities under two assumptions.
First, to automatically detect discontinuities we need an architecture that can exploit the correlations between the values of the solution for all spatial locations x 1 , . . . , x N . Second, the architecture has to be flexible enough to capture different behaviors of the solution at different regions of the domain. To this end, we propose the use of encoder-decoder GRUs 2 for predicting the solution at all locations at once, with the use of a recent machine learning tool known as attention mechanisms 52 .
Our approach presents several advantages compared to traditional simulators: (i) Instead of using just neighboring cells' information to calculate u as in numerical methods, our architecture uses the complete encoded sequence input of the grid to obtain u i , allowing us to capture non-local relationships that numerical methods struggle to identify. (ii) The computer time for the forward pass of neural networks models is linear with respect to the number of cells in our grid. In other words, our method is a faster alternative with respect to traditional numerical methods of solving PDEs such as finite difference. Figure 1 shows an outline of the proposed architecture. We start feeding the input pair (t, M) to a single fully connected layer. Thus, we obtain h 0 the initial hidden state of a sequence of N GRU blocks (yellow). Each of them corresponds to a spatial coordinate x i which is combined with the previous hidden state h i−1 inside the block.
This generates a set of vectors y 1 , . . . , y N which can be understood as a representation of the input in a latent space. The definitive solution u (we omit the subindex θ for simplicity) is reached after a new sequence of GRU blocks (blue) whose initial hidden state d 0 is initialized as h N to preserve the memory of the system.
In addition to the hidden state d i , the i-th block g i is fed with a concatenation of the solution at the previous location and a context vector, that is How the context vector is obtained is one of the key aspects of our architecture, since it will provide the PINN with enough flexibility to fit to the different behaviors of the solution depending on the region of the domain. Inspired by 51 , we introduce an attention mechanism between both GRU block sequences. Our attention mechanism is a single fully connected layer, a, that learns the relationship between each component of y j and the hidden states of the (blue) GRU sequence, Then, the rows of matrix E are normalized using a softmax function as (10)  The coefficients α i,j can be understood as the degree of influence of the component y j in the output u i . This is one of the main innovations of our work to solve hyperbolic equations with discontinuities. The attention mechanism automatically determines the most relevant encoded information of the full sequence of the input data to predict the u i . In other words, attention mechanism is a new method that allows one to determine the location of the shock automatically and provide more accurate behavior of the PIANN model around this location. This new methodology breaks the limitations explored by other authors 39,45,46 since is able to capture the discontinuity without specific prior information or the regularization term of the residual. This is the first paper to use attention mechanisms to solve PDEs. Furthermore, we show that attention mechanisms are particularly effective in capturing the challenging nonlinear features of hyperbolic PDEs with discontinuous solution.

Results
To test the proposed PIANN methodology, we perform a set of numerical experiments where we train our network to learn the solutions of the Buckley-Leverett problem (1)-(3) for a wide range mobility ratios. We compare our numerical predictions with the analytical solution given in Eq. (5).
The training set is given by a grid  Figure 2 shows the residual value for the testing dataset for the different epoch for M = 4.5 . We can observe a fast convergence of the method and a cumulative value of the residual smaller than 10 −4 after a few epochs. This demonstrates that we are minimizing the residuals in Eq. (1) and subsequently solving the the equation that governs BL. Figure 3 shows the comparison between the analytical solution (red) and the solution obtained by our PIANN (blue) for different M used during training. Top, middle and bottom rows correspond to M = 2 , M = 48 and M = 98 , respectively, and the columns from left to right, correspond to different time steps t = 0.04 , t = 0.20 , and t = 0.40 , respectively. We can distinguish three regions of interest in the solution. Following increasing x coordinates, the first region on the left is one where the water saturation varies smoothly following the rarefaction part of the solution. The second region is a sharp saturation change that corresponds to the shock in the solution and the third region is ahead of the shock, with undisturbed water saturation values that are still at zero. For all cases, we observe that the PIANN properly learns the correct rarefaction behavior of the first region and approximates the analytical solution extremely well. In the third region, the PIANN also fits to the analytical www.nature.com/scientificreports/ solution perfectly and displays an undisturbed water saturation at zero. As for any classical numerical methods, the shock region is the most challenging to resolve. Around the shock front, the PIANN seems unable to perfectly propagate a strict discontinuity, and the represented behavior is a smeared shock that is smoothed over and displays small non-monotonic artifacts upstream and downstream of the front. The location of the shock is, however, well captured, suggesting that the network has correctly learnt the conservative nature of the BL equation. The non-monotonic behavior is reminiscent of the behavior observed in higher-order finite difference or finite volume methods, where slope-limiters are often used to correct for the non-physical oscillations.
The PIANN needs to learn where the shock is located in order to fit differently to both sides of it. This is the role of the attention mechanism of our architecture. On top of each figure we have visualized the attention map introduced by 51 for every timestep. These maps visualize the attention weight α i,j to predict the variable u i . We observe that in all cases the attention mechanism identifies the discontinuity, water front, and subsequently modifies the behavior of the network in the three different regions described above. This shows that attention mechanism provides all the necessary information to capture discontinuities automatically without the necessity of training data or a prior knowledge. Finally, it is important to note that attention weights of the attention mechanism are constant when the shock/discontinuity disappears.
We test the behavior of our methodology to provide solutions for BL at points within the training range of M: M = 4.5 and M = 71 . In other words, we want to check the capability of our PIANN model to interpolate solutions. Figure 4 shows that our PIANNs provide solutions that correctly detect the shock. Figure 5 shows the absolute error for both mobility ratios. As we can observe, the absolute error is zero after the shock and the average absolute error is smaller than 10 −3 .
In addition to interpolation, we have studied the behavior of PIANN for extrapolation, in other words, to solve BL for mobility ratios outside the initial training set, M ∈ [2, 100] . Specifically, we explore the network predictions for values of M outside the training set: M = 140 , M = 250 and M = 500 . Figure 6 compares the predicted and the exact solutions of the PDEs. As expected, the accuracy degrades when M is far from the original training set. We observe that our method predicts the behavior of the shock and accurate solution for M = 140 and M = 250 . However, the shock is totally missed for M = 500 and, as such, retraining the model is recommended for larger values of M. It is important to note that the results show that our method is stable and converges for the different cases.
We test how the neural residual error progresses based on different t and x resolutions. Results are shown in Table 1, and demonstrate that our PIANN obtains smaller residuals when the resolution of the training set increases. However, we observe that changes in the residual are not highly significant. This is an advantage with respect to traditional numerical methods in computational fluid dynamics, where smaller values of t are necessary to capture the shock and guarantee convergence and stability.
Finally, we have compared the results with central and upwind finite difference schemes for the term of the vector of fluxes. The first-order upwind difference introduces a dissipation error when applied to the residual of the Buckley-Leverett equation, which is equivalent to regularizing the problem via artificial diffusion. Figure 7 shows that both approaches present similar results respect to the analytical solution. The fact that both central and upwind differences yield similar predictions is important, because it suggests that the proposed PIANN approach does not rely on artificial dissipation for shock capturing.

Discussion
In this work, we have introduced a new method to solve hyperbolic PDEs. We propose a new perspective by focusing on network architectures rather than on residual regularization. We call our new architecture a physics informed attention neural network (PIANN). www.nature.com/scientificreports/ PIANN's novel architecture is based on two assumptions. First, correlations between values of the solution at all the spatial locations must be exploited, and second, the architecture has to be flexible enough to identify the shock and capture different behaviors of the solution at different regions of the domain. We have proposed an encoder-decoder GRU-based network to use the most relevant information of the fully encoded information, combined with the use of an attention mechanism. The attention mechanism is responsible for identifying the shock location and adapting the behavior of the PIANN model.
The loss function of PIANNs is based solely on the residuals of the PDE, and the initial and boundary conditions are introduced as hard constraints in the architecture. As a result, PIANN's training aims only at minimizing the residual of the PDE; no hyperparameters are needed to control the effect of initial and boundary conditions on the solution.
We have applied the proposed methodology to the non-concave flux Buckley-Leverett problem, which has hitherto been an open problem for PINNs. The experimental results support the validity of the proposed methodology and conclude that: (i) during training, the residuals of the equation decrease quickly to values smaller than 10 −4 , which means that our methodology is indeed solving the differential equation, (ii) the attention mechanism automatically detects shock waves of the solution and allows the PIANN to fit to the different behaviors of the analytical solution, and (iii) the PIANN is able to interpolate solutions for values of the mobility ratio M inside the range of training set, as well as to extrapolate when the value of M is outside the range. However, we observe that if M is too far away from range of the training set, the quality of the solution decreases. In that case, retraining of the network is recommended. (iv) We observe that the residuals decrease when the resolution of the training set increases. However, the change in the residuals is not highly significant. This is advantageous  www.nature.com/scientificreports/ with respect to traditional numerical methods where small values of t are needed to capture the shock and guarantee convergence and stability.

Conclusion
In conclusion, the proposed methodology transcends the current limitations of deep learning for solving hyperbolic PDEs with shock waves, and opens the door to applying these techniques to real-world problems, such as challenging reservoir simulations or carbon sequestration. However, the scaling up the proposed method for 3D problems will require to compute attention mechanism very intensive for memory, attention mechanism scale quadratically. A possible approach to overcoming this problem is to use a sparse representation of neural network, use attention mechanism to determine the most relevant time sequences or use simpler attention mechanism that has a linear dependency respect to number of cells. Finally, it is plausible that this method could be applied to model many processes in other domains which are described by non-linear PDEs with shock waves. The extent to which PIANNs could be applied in the more general space of differential equations is left to future research.  www.nature.com/scientificreports/